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ABSTRACT 

We perform a suite of cosmological simulations in the ACDM paradigm of the formation of the 
first structures in the Universe prior to astrophysical reheating and reionization (15 ^ z < 200). 
These are the first simulations initialized in a manner that self consistently accounts for the impact 
of pressure on the rate of growth of modes, temperature fluctuations in the gas, and the dark matter- 
baryon supersonic velocity difference. Even with these improvements, these are still difficult times 
to simulate accurately as the Jeans length of the cold intergalactic gas must be resolved while also 
capturing a representative sample of the Universe. We explore the box size and resolution requirements 
to meet these competing objectives. 

Our simulations support the finding of recent studies that the dark matter-baryon velocity difference 
has a surprisingly large impact on the accretion of gas onto the first star-forming minihalos (which 
have masses of - 10^ Mq). In fact, the halo gas is often significantly downwind of such halos and 
with lower densities in the simulations in which the baryons have a bulk flow with respect to the 
dark matter, modulating the formation of the first stars by the local value of this velocity difference. 
We also show that dynamical friction plays an important role in the nonlinear evolution of the dark 
matter-baryon differential velocity, acting to erase this velocity difference quickly in overdense gas as 
well as sourcing visually- apparent bow shocks and Mach cones throughout the Universe. 

We use simulations with both the GADGET and Enzo cosmological codes to test the robustness 
of these conclusions. The comparison of these codes' simulations also provides a relatively controlled 
test of these codes themselves, allowing us to quantify some of the tradeoffs between the algorithms. 
For example, we find that particle coupling in GADGET between the gas and dark matter particles 
results in spurious growth that mimics nonlinear growth in the matter power spectrum for standard 
initial setups. This coupling is alleviated by using adaptive gravitational softening for the gas. In a 
companion paper, we use the simulations presented here to make detailed estimates for the impact 
of the dark matter-baryon velocity differential on redshifted 21cm radiation. The initial conditions 
generator used in this study CICsASS can be publicly downloaded. 

Subject headings: cosmology: theory — first stars — galaxies: high redshift - stars: Population III - 
galaxies: formation 



1. INTRODUCTION 

For the first hundred million years after the Big Bang, 
the gas distribution in the Universe as well as its ther- 
mal state can be accurately calculated by solving a set 
of linear differential equations. However, between the 
redshift s of 100 and 10, the inhomogeneties in much of 
the cosmic gas went nonlinear. Eventually, deep enough 
potential wells for the primordial gas to cool and form 
stars developed, and the Universe transitioned to a vastly 
more complex system in which stars abound and their 
radiative backgrounds impacted all of the baryonic mat- 
ter. In principle, it is possible to understand perfectly 
the evolution of gas and dark matter before stellar radi- 
ation backgrounds impacted all matter, at z > 20, us- 
ing a combination of linear theory and (once nonlinear 
structure begins to form) numerical simulations. Never- 
theless, important questions about the evolution of the 
Universe during this virgin epoch remain unanswered. 

For example, it is unclear whether weak structure for- 
mation shocks would have significantly heated the cos- 
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mic gas (Gnedin &: Shaver| 2004 [Furl anetto fc Loeb 

CK! 



[2004 ) . In the absence of shocks, the intergalactic medium 
(TGM) is anticipated to have been kinetically cold prior 
to reheating by astrophysical sources, reaching 10 K at 
z = 20, and with the temperature cooling adiabatically 



as (1 



However, even 0.3 km s flows would have 



been supersonic for a gas temperature of 10 K, and su- 
personic motion is likely to source shocks and entropy 
generation. 



In addition, recently Tseliakhovich & Hirata (2010) 
demonstrated that in most places in the early Universe 
the baryons and dark matter were moving supersonically 
with respect to one another. At the time of recombina- 
tion, the cosmic gas was moving with respect to the dark 
matter at an RMS velocity of 30 km s~^ and in a coher- 
ent manner on < 10 comoving Mpc separations. These 
initial velocity differences translate into the dark matter 
moving through the gas with an average Mach number 
of A^bc ^ 1-7 over 15 < z < 150, but with a standard de- 
viation in Mach number between different regions in the 
Universe of 0.7. As with structure formation, such su- 
personic motion may source shocks, generating entropy 
and reheating the universe. 

In this paper, we simulate the evolution of matter in 
the Universe prior to when radiation backgrounds gener- 
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ated by stars became important sources of heating. This 
time has been the focus of many prior studies of th e 
first stars (e.g., [Abel et al.||2QQ2[ [Bromm et al.|[2QQ2| . 
We discuss the simulation box size and resolution re- 
quirements to simulate these times accurately, and we 
add several improvements to standard methods for ini- 
tializing cosmological simulations so that our simulations 
are initialized with full linear solutions for the growth of 
structure. For example, in contrast to prior studies, our 
initial conditions self-consistently account for the impact 
of gas pressure on the growth of modes as well as include 
fluctuations in the gas temperature (an improve me nt em- 
phasized as impor tant in ,Naoz fc Barkana|2QQ5 and Naoz 
fc BarkanaJ ( [20071 )). 



Yue 0.24, a nd 



A tew prior studies have investigated the impact of the 
supersonic motion of the gas relative to the dark mat- 
ter on the formation of the firs t gas-r ich hal os and on 



the first star s (Maio et al . 2011 Stacy et al. 2011 Greif 



et al. 



2011b 



et al 



Interestingly, some ot 

these studies hnd this motion lias a dramatic impact on 
the formation of the first stars. However, the relative 
velocity in all of these simulations was incorporated by 
boosting the velocity of the gas at the onset of the sim- 
ulation, which we show misses much of the impact of 
this supersonic motion on the linear growth of structure. 
The simulations in this study are the first to use a con- 
sistent linear theory to initialize these differential flows. 
Our simulations enable us to more rigorously test these 
claims as well as to investigate other manifestations of 
such cosmic flows. 

The Universe during the 'Dark Ages' - times before 
stars reionized and reheated the Universe - is observable 
via the redshifted 21cm line in absorption against the 
cosmic microwave background. Several collaborations 
are currently developing instruments to detect this era 
(LEDA, DARE, and LOFAR^). The strength of the sig- 
nal, especially on large scales, is intimately related the 
both ther mal history of the gas as well as the star forma- 
tion rate (padau et al.|1997i|Furlan etto'2006>'Furlanetto" 
et al. |200D|. In a companion paper (McQuinn fc OTeary 



2Ul2; hereafter Paper II), we will address the observa- 
tional signatures of this era in the redshifted 21cm ab- 
sorption signal, specifically focusing on the impact of the 
relative velocity between the baryons and dark matter. 

This paper is organized as follows: Section |2] elucidates 
the characteristic scales and physical processes that 
affect the evolution of intergalactic gas at 15 < z < 200. 
Section |3] discusses considerations relevant to simulating 
these early cosmic times as well as the details of our 
initial conditions generator. This section also compares 
cosmological simulations of th e early U niverse run using 
both the GADGET ( Springel et al .l 2001) and Enzo 
|0'Shea et al.||2004[ ) codes. Section [s] describes the 
important roll of dynamical friction in the non-linear 
evolution of structure formation with a baryonic stream- 
ing velocity, v\yc- Section [6] uses these simulations to 
characterize the properties of this cosmic epoch as 
well as the impact of v\)c on structure formation. This 
study assumes a flat ACDM cosmological model with 
= 0.27, = 0.73, h = 0.71, ag = 0.8, = 0.96, 
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0.046, consistent with recent 
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measurements (Larson et al. [2011 ). 
z = dark matter traction as lie = ^ 
subsequently abbreviate proper Mpc as pMpc, and Mpc 
and kpc will be reserved for comoving lengths. Some of 
our calculations use the Sheth-Tormen mass function, 
for which we adopt the paramete rs p = 0.3, a = 0.75, 
A = 0.322 dSheth fc Torinen|[2QQ2l ) . 



2. CHARACTERISTIC SCALES IN THE 
POST-RECOMBINATION AND PRE-REIONIZATION 
UNIVERSE 

The focus of this paper and Paper II is on the era after 
the gas thermally decoupled from the CMB and before 
it was reheated by astrophysical sources. This period is 
anticipated to have occurred over 15 < z < 200 ( Peebles] 
pl93, Gnedin fc Shaver|2004[|Furlanetto|2006^ . It should 
be possible to model most aspects of these pristine times 
from the cosmological initial conditions alone. The tem- 
perature of the gas at the cosmic mean density cooled 
adiabatically with the expansion of the Universe during 
this period, decreasing with redshift as (1+z)^ and reach- 
ing a temperature of 10 K at z = 20. In fact, the vast 
majority of the gas likely existed near a single adiabat: 
It is unclear when or if structure formation shocks would 
have contributed significantly to the entropy of the IGM 
(a question investigated here), and only in the most mas- 
sive, rarest halos was the gas able to cool by molecular 
or atomic transitions. In addition, only at the end of 
this adiabatic period was more than a percent of the gas 
bound to dark matter halos. In particular, the fraction 
of matter that had collapsed into dark matter halos with 
masses > 3 x 10^ Mq - those massive enough to over- 
come pressure and retain gas fNaoz et al. 2009) - at 
z = 15, 20, and 30 was 0.03, 0.005, and 6 x lO^Taccord- 
ing to the Sheth-Tormen mass function. These numbers 
become 0.008, 7 x lO""^, and 1 x 10"^ for > lO^M© halos 
- halos massive enough to cool via molecular h ydrogen 
transitions and host stars ( Tegmark et al.|p-997 ). 

Figure [T] shows the scales most relevant to the said 
epoch in comoving coordinates: the virial radius of a 
10^ Mc^ halo and of one with 10^ A f^^, the Jeans' and 



Filtering length as defined in Naoz fc Barkana (2 007| ), as 
well as the radius of a sphere in which the RMS linear 
density contrast, cr^, equals 0.1, 0.2, 0.3, 0.4, and 0.5. 
The virial radii of halos that can retain their gas and 
form stars, the Jeans' length, and the Filtering length 
are all ~ 1 — 10 kpc. In fact, the comoving Jeans' length, 
is nearly constant over 15 < z < 200, and equal to 
5 kpc with the scaling Rj oc (1 -h^)^/^(l + ^)^/^, where 6 
is the matter overdensity. The comoving Filtering length, 
which is the analogue of the Jeans ' lengi:th at the mean 
density of an expandin g universe ( Gnedin fc Hui|p"998 ' 



Naoz fc Barkana 2007), also has a weak dependence on 
redshift.^ The vertical line segments in Figure [l] show 



^ Interestingly, Na oz et al.| (^009 ) found that the linear-theory 
Filtering mass characterizes tiie minimum mass halo that can over- 
come pressure and retain its gas; smaller mass halos are largely 
devoid of gas. We think there is a simple explanation for why 
the Filtering mass sets the minimum mass of a halo that can re- 
tain its gas. For adiabatic collapse (appropriate for almost all gas 
at the specified epoch), the Jeans' mass increases with density as 



In addition, for gas in the Hubble flow, the characteristic 
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Fig. 1. — Characteristic comoving scales in the Universe prior to 
astrophysical reheating: the virial radius of 10^ Mq and 10^ Mq 
halos, the Jeans' length for mean density gas, the Filtering length, 
and the radius of a sphere with RMS linear density fluctuation 
of the specifled value. Each vertical bar represents one of the sim- 
ulations employed in this study, with its height stretching from the 
mean interparticle distance (or base-grid mesh width) to the sim- 
ulation's box size. The leftmost bar represents a simulation with 
{0.1/ h Mpc, 256^ gaseous resolution elements}, and the rightmost 
one with {1/h Mpc, 768^ gaseous resolution elements}. 

the range of scales that are captured by the simulations 
we run for this study and discussed in ^ 

3. SIMULATING THE DARK AGES 

3.1. Simulation Parameters and Initial Conditions 

There are a few hurdles that must be overcome in 
order to accurately simulate the formation of the first 
structures in the Universe and the impact of the dark 

scale above which gas fragments is the Filtering mass. For realis- 
tic thermal histories, the Filtering mass is even smaller than the 
Jeans' mass, at least at the time when a region that virializes at 
z ~ 20 decoupled from the Hubble flow. Therefore, collapsing gas 
formed its smallest unit when it was in the Hubble flow, and these 
units were unable to fragment further as they collapsed and in- 
stead were dragged into the dark matter potential wells as they 
formed. Following the same logic, the flltering mass should be a 
poor approximation for the characteristic halo mass that is able 
to retain its gas after reionization. In this case, the gas was pho- 
toheated to 10^ K by reionization prior to collapse. This initial 
temperature floor results in the collapse being non-adiabatic and 
instead being better approximated as isothermal once the gas had 
collapsed to moderate overdensities such that it has adiabatically 
heated to ~ 2 — 3 x 10^ K (temperatures where collisional cooling is 
extremely efficient). For isothermal gas, the Jeans' mass scales as 
p~^l^\ collapsing regions fragment into smaller and smaller clumps 
as their density increases. In support of the hypothesis, Hoeft et al!] 
|2006 ) and Okamoto et al. (2008 ) find using numerical simulations 
that this characteristic halo mass that retains gas after reionization 
is generally substantially smaller than the Filtering mass. In addi- 
tion, this also provides an explanation for why studies using cosmo- 
logical simulations flnd that the sizes of overdense gaseous clumps 
in the photoionized IGM are approximately set by the Jeans' scale 
for the clump's current density rather than the Jeans' scale at the 
cosmic mea n density (|Schaye||200lj iMcQuinn et al.||2011| lAltayl 
|et al.|2011| >. 
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Fig. 2. — Predicted halo mass function, dn/drrih, for simulations 
with box sizes that span the range studied here (as well as a hy- 
pothetical 2 Mpc/h box). Shown is the expected mass function, 
which we set to zero when mhdn/dmhL^^^ < 1. The rightward 
pointing arrow indicates the halo masses at which inter-halo gas 
has the potential to cool via molecular emissions and form stars 
(using the criterion -Ucir > 3.7 km s"-*^; e.g., (Fialkov et al.|201l] ). 



matter-baryon velocity differential: 

First, a simulation of these times needs to capture the 
pressure smoothing scale of the gas (or ~ 5 comoving 
kpc between 20 < z < 200; Fig. [T]) while also being 
large enough so that the box-scale modes are still lin- 
ear. Otherwise, it does not properly capture the gas 
physics and/or is not representative of the Universe. As 
illustrated in Fig. [l] capturing a representative volume 
becomes increasingly difficult with decreasing redshift 
around z ^ 20 because the nonlinear scale rapidly moves 
to larger scales with decreasing redshift owing to the 
near scale-invariance of Mpc-scale density fluctuations. 
A non-representative box size has a dramatic effect on 
the halo mass function (Barkana & Loeb 2004). Figure [2] 
shows the range of halo masses that different simulation 
box sizes capture.^ A 1 Mpc/h box is needed to sta- 
tistically capture the z = 20 halo mass function at the 
factor of < 2-level at the mass threshold that can cool 
by molecular hydrogen transitions, and an even larger 
box size is required to meet this requirement to higher 
redshifts. 

Second, the effects of the baryonic streaming velocity 
on the growth of Jeans '-scale perturbations is important 
over a broad range of redshifts after the baryons ther- 
mally decoupled from the CMB ( Tseliakhovich & Hirata 



mally decoupled trom tne LjMd ( 1 seliaknovicn 6^ Mirata 
|201 0|. Fi g ure j 3| shows solutions to the linearized equa- 
tions ( Al - |A4| m Appendix [A|) using initializations of ^'bc 
common m the numerical literature as well as the full lin- 

^ The halo mass function as a function of box size in Fig- 
ure |2] is calculated by multiplying the Sheth-Tormen mass function 

by npsinihl^a'^^ - crf^^^^2)/nps(mh\am), where nps(m^,crx) 

is the Press- Schechter mass function at mass rrih given ax, the 
RMS density con trast in a sphere of rad ius X. This prescription 
was motivated in [Barkana Sz Loeb| ( |2004| . 
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Fig. 3. — Investigation of different approximations for initializing 
numerical simulations with ^'bc• Each curve is the linear growth 
factor at z = 30 of the dark matter (top panel) and the gas (bottom 
panel), initialized so that it equals unity at z = 1000, and for 
modes with -UbcCOS^fc = 3 km s"-*^ at z = 100 (A4bcCOS(/)fc = 1.8) 
unless labelled otherwise. All previous simulations added a relative 
dark matter-baryon velocity on top of the standard cosmological 
initial conditions at the initialization redshift of the simulation, 
either z = 100 or 2; = 200. As can be gauged from comparing 
to the sold curve, which represen t the full solution as discussed 
in [Tseliakhovich Hirat aj (12010^, to the long and short dashed 
curves, such an mitialization misses much of the effect of vi^c in 
linear theory. In addition, previous simulations had not included 
pressure self-consistently in the initial conditions. If the baryons 
are initialized at 2; = 100 with A4hc = 1-8 and such that they 
trace the dark matter distribution as in most previous studies, for 
cos = 1 this leads to the cyan dot-dashed curve in the top panel. 



ear theory initialization used in this paper. Each curve 
is the hnear growth factor at z = 30 of the dark mat- 
ter (top panel) and the gas (bottom panel), initialized so 
that it equals unity at z = 1000. As can be gauged from 
comparing to the sold curve (which represent the initial 
conditions used in this study where vi^c is consistently 
followed) to the long and short dashed curves (which ap- 
proximate the initial conditions use d in Maio et al.|2011[ 
[Stacy et al.(2QTT||Greif et al.|2011b[ and |Naoz et al"|2011 
where v\yc is incorporated only at the onset of the simu- 
lations), simply boosting the velocity of the baryons at 
z = 100 or 200 misses much of the linear effect on 5^ 
and 6c. The dotted curve is the case with '^bc = 0. In 



addition, the vast majority of studies have assumed that 
the baryons also trace the dark matter, which results in 
the baryons streaming out of the potential wells of the 
dark matter early on in the simulations and leads to the 
dot-dashed curve in the top panel. 

Lastly, simulations exploring the high-redshift Uni- 
verse cannot be initialized at similar redshifts to those 
that are used to understand lower redshifts. A halo that 
collapses at z = 25 has a linear over dens ity of 0.44 at 
z — 1 00 in the spherical collapse model (Gunn & Gott 
1972) and so the error in the overdensity from using 



only Eularian (Lagrangian) linear theory is uncomfort- 
ably large, 52% (11%), with this error resulting in the 
str uctures being less bou nd and collapsing at later times. 
See Crocce et al. ( |2006 | for more quantitative determi- 
nations of this error. Initial conditions that use 2"^ order 
perturbation theory would reduce this error, but higher 
order solutions that account for gas pressure have only 
been developed for toy cases (Shoji & Komatsu 2009 ). At 
the same time, the particle noise in simulations initialized 
at too high of a redshift dominates over the cosmologi- 
cal clustering. Thus, a balance must be reached between 
initializing at a redshift where 1®* order Lagrangian per- 
turbation theory is accurate and avoiding particle noise. 
Prior studies of y^c (a - s well as mo st studies of the first 
stars; JAbel et al ||2Q02| [Greif et al.p Ollb) typically used 
Zi ^ loo and l""^ order Lagrangian perturbation theory. 
We favor higher redshifts in this study, with Zi ~ 200 or 
400, since this choice does not significantly increase the 
computational requirements of the simulations nor add 
much extr a noise to the power spectrum of the density 
field (p^. 

To initialize our simulations, w e adopt the approxima- 



tion in Tseliakhovich fc Hira ta 
tivated in 5.3.1 in |Hu||l995 



(. 2010) (and further mo- 
dlEisenstein & Hu||l998] 



that the baryons dec ouple^n stantaneously at z = lOQi 
and solve equations ( A1|A4 ) for the growth of matter 
fluctuations in linear theory to redshift Zi. These 1^* or- 
der differential equations are initialized with the CAMB 
Boltzmann code transfer function^ and its time deriva- 
tive. We find that this approximation reproduces the 
evolution of the matter power in CAMB excellently for 
^bc = 0- To relate the linear theory solution to particle 
displacements, we use that the linearized displacement 
^ from Lagrangian position q is related to the Eular- 
ian linear perturbation theory overdensity, Se, via the 
relationship V • ^ = —5e for an irrotational flow (th e 



x = q^^ = q - V^, 

v = —a V^, 



"Zel'dovich approximation"; e.g., Padmanabhan 1993). 
It follows that the linear displacement and its velocity 
are given by 

(1) 

(2) 

where V^^ = Se- These equations set the displace- 
ments that are used from the initial positions set by a 
glass file (using cloud-in-cell interpolation from a grid 
with each dimension equal to the cube root of the par- 
ticle number).^ For our SPH simulations, we typically 

^ http://camb.info/ 

^ Linear order Lagrangian perturbation theory (the Zel'dovich 
approximation) smooths out small scales, resulting in the dark mat- 
ter power spectrum at lowest order being given by 



Pz{k) = PEik)exp[-ajik^/2], 



(3) 



where Pe is the linear theory Eularian density power and aj^ = 
(37r)-2 IjI""^^^ dk Psik). Thus, this theory wih be valid at /c < cr~^ 

such that Pz(k) = Pe- In a 1 Mpc box, (Tr = 0.3 G{z) Mpc, where 
G{z) is the growth factor. To illustrate this difficulty, for a 0.1 Mpc 
box, fc* = (7^ = 3400 Mpc"-*^ at 2; = 200, which is comparable 
to the Nyquist frequency, kjy = 8000 (A^/256) Mpc~^ Here, 
is the cube root of the particle number. For a 0.5 Mpc box, fc* 
becomes fc* = 900 Mpc~^ and for a 1 Mpc box, fc* = 500 Mpc~^ 
These numbers may be prohibitive if we were attempting to capture 
the smallest scales (and smallest dark matter halos) in our box. 
However, fc^r is a scale that is generally buried in the particle noise. 
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also use a differe nt glass file to initial ize the baryons as 
recommended in Yoshida et al. (2003) to avoid spurious 
couplings between ditterent particle types, although we 
found this resulted in increased coupling in some simu- 
lations (see Appendix B). 

The major improvements of our algorithm over other 
cosmological initial conditions generators are: 

• The baryons and dark matter have distributions 
and velocities that are consistent with linear theory, 
including the impact of pressure. Other commonly 
used initial conditions generators assume that the 
two components trace each other (the GADGET 
publicly available initial conditions code), that the 
dark matter and baryon velocities are proportional 
to each components' respective transfer function 
(the Enzo distribution "inits" code), or that the 

two compo nents have the same velocity ( Yoshida^ 

et al. |2003[ ). Other codes assume some variant ot 
= Qrn{z)^^^H which is not valid at scales 
where pressure is important nor when radiation im- 
pacts cosmological expansion. 



• We include the effects of radiation in the initial den- 
sities and velocities. Radiation is also included in 
the background evolution in the simulation them- 
selves. For both, we assume that the three species 
of neutrino are relativistic at all redshifts. Radi- 
ation impacts the rate of growth at the 10% level 
for simulations starting at z ~ 300. 

• The mean temperature and electron density are ini- 
tialized with the values calculat ed with the REC - 
FAST recombination code (Seag er et al. 1999 ).^ 
Our initial conditions also include fluctuations in 
the gas temperature as ca lculated from linear the- 
ory, an improvement that Naoz & Barkana (2007) 
stressed as important.^ 



• When a relative bulk velocity between the dark 
matter and baryons ('^bc > 0) is required, we use 
linear solutions that self-consistently incorporate 
v\)c- This improvement is emphasized in the en- 
suing discussions. 

This initial conditions generator, the Cosmological Ini- 
tial Conditions for AMR and SPH Simulations (CIC- 
sASS) can be downloaded at astro . berkeley . edu/$\| 
rsimSmmcquinn/ codes , 

The criterion that needs to be satisfied is /c* ^ kp, where the 
Filtering wavevector kp is ~ 500 Mpc"-*^. In all of the simulations 
employed here, this condition is at least weakly satisfied. Zoom in 
simulations of the first stars, where fcmin can be much smaller than 
in our simulations, should be wary of this deficiency of 1^^ order 
Lagrangian perturbation theory. 

^ Compton cooling is also self-consistently included in the simu- 
lations. GADGET and Enzo first stars calculations appear to use 
Case A recombination coefficients, but Case B (and calibrated for 
low temperatures) is a better choice and is required to obtain the 
correct evolution in the electron fraction and, hence, the correct 
thermal history. 

^ The prefactor for converting temperature units in the Enzo 
code appears to have a typo (1.88e6 should read 1.8182e6), result- 
ing in an absolute error in the temperature of ~ 3.4%. We have 
not corrected our calculations for this error. 



3.2. Numerical Simulations 

We use these initial conditions with th e cosmologi- 
cal cod es GADGET3 (Springel et al. 2001) and Enzo 
v2.1.1 (O'Shea et al. 2004) with the default piecewise 
parabolic method to solve the equations of hydrodynam- 
ics. GADGET solves the equations of fluid dynamics 
with the smooth particle hydrodynamics (SPH) method, 
whereas Enzo is a grid code with adaptive mesh refine- 
ment (AMR). Furthermore, for gravity GADGET uses a 
the tree-particle mesh grid whereas Enzo a nested parti- 
cle mesh grid. Both the GADGET and Enzo codes have 
been shown to conserve entropy at the part in lO QO-level 
for th e expansion of a homogeneous Universe (O'Shea 
Fna1f2005). Such accuracy is unusual among hydro- 
dynamics codes, and owes to the entropy-conserving for- 
malism of GADGET and the 3^^-order accurate in space, 
2"^-order in time Riemann solver employed by Enzo. 
In addition, we require the codes to capture the weak 
shocks that develop from the mildly supersonic flows in 
the simulations. Shock capturing is a strength of the 
Enzo algorithm but a potential weakness of GADGET 
(and SPH codes in general), which does not explicitly 
capture shocks. Hence, our work includes a direct com- 
parison between these two codes. 

The GADGET simulations that were run include {box 
size in Mpc//i, gas particle number} of {0.1, 128^}, 
{0.1, 256^}, {0.2, 256^}, {0.2, 512^}, {0.5, 512^}, and 
{1, 768^}, ah initialized at a redshift of Zi = 200.^^ Our 
GADGET simulations were run using adaptive gravita- 
tional softening of the gas particles to reduce the amount 
of gas-dark matter particle coupling (Appendix [B|). At 
each box size and particle number, baryon streammg ve- 
locities, 'Ubc, that result in Mach numbers of Mhc = 
and 1.9 during the Dark Ages were both seperately run. 
(Note that, at the scales captured in our simulation 
boxes, the baryons were initially moving coherently with 
velocity vi^d^i) as a uniform wind. The spatial distribu- 
tion of 'Ubc is determined by the photon-baryon acoustic 
physics for which the spatial fluctuations are damped on 
scales captured by our simulation boxes by Silk Damp- 
ing.) These choices allowed us to explore the sensitivity 
of our results to resolution and sample variance. Note 
that all simulations resolve to varying degrees the Jeans' 
scale (Fig.jl]), and the dark matter (gas) particle mass in 
the {0.5, 512^} simulation is 82 (17 M©). We also 
reran a sample of these simulations with (1) A^bc = 0-6, 
(2) Aihc = 3.8, (3) a more sophisticated separate ar- 
tificial viscosity implementation, (4) with a fixed gravi- 
tational softening length for the gas, and (5) Zi = 400. 
Finally, we do not adopt a standard practice in early 
Universe simulations of increasing os to compensate for 
missing large-scale power. We found this artifice makes 
the results difficult to interpret. 

The Enzo simulations we ran include runs with a uni- 
form grid (with no AMR) with {box size in Mpc//i, grid 
size} of {0.1, 256^}, {0.2, 256^}, {0.2, 512^}, and 
{0.5, 512^} with both Mhc = and 1.9 (and in a couple 
cases 3.8). We also ran AMR runs with {0.2, 256^} and 
{0.5, 512^} with 4 levels of adaptivity for both the hy- 
dro grid and for the gravity grid. We refined on gas den- 

The 768^ simulations use the more sophisticated viscosity 
implementation of ^Morris &: Monaghan^ ( ^1997j . 



Fig. 4. — Slices of log(l + S^) through z = 20 snapshots of the {0.2 Mpc//i, 2 x 256^ resolution element} GADGET (top panels) and 
Enzo (bottom panels) simulations. The top left, middle, and right panels correspond respectively to the cases of -Ubc = [A^bc = 0]: 
-Ubc = 30(2/1000) km s~-^ [J^hc = 1-9], and t'bc = 60(2/1000) km s"-*^ [A^bc = 3.8]. Dark regions represent overdensities and light 
underdensities, and the contrast is the same in all of the panels. These simulations were initialized with the same random numbers and 
such that the baryons were flowing to the right in the cases with nonzero ^'bc• The Enzo simulations used a flxed uni-grid with 512"^ cells. 

lations and all of the adaptive runs of the Enzo simu- 
lations, molecular hydrogen formation was not included. 
Molecular hydrogen is the only active coolant at the gas 
temperatures in our simulations, and its absence prevents 
the runaway collapse into stars of dense gas in the first 
halos. 

We have subjected our simulations to a battery of tests 
in order to confirm the robustness of our results. For both 
GADGET and Enzo, we looked at the grid size, box size, 
maximum time-step size, number of particles/cells, and 
frame of reference for the relative velocity of the baryons 
and dark matter on the grid. We found that the statistics 
we considered, such as the volume- averaged temperature 
and the power spectrum of the matter, were robust to 
these choices aside from grid size and box size (and we 
show these dependencies in ensuing plots). For Enzo, we 
also investigated different Riemann solvers (and saw no 
differences) and the number of AMR levels. For GAD- 
GET, we looked at the impact of using an adaptive ver- 
sus a fixed gravitational softening length for the g as, as 
well as different ways of staggering the initial distribution 
of particles and different artificial viscosity implementa- 
FiG. 5.- Slices through the GADGET (top panel) and Enzo (b^^^^ tions.^^ We found that, when using an adaptive gravita- 
tom panel) simulations. Shown is the same slice through log(l + Obj _ , • 1 n 7 

in the left panels in Fig. [3 (the case with TWbc = 3.8), except that Clonal softenmg length, there was excellent convergence 

only the top half of Fig.^ is shown and the contrast is linear for 
0-5 ^ ^ 2.0. The dashea, rotated 'V in the bottom panel shows 
the opening angle of a bow shock for A^bc = 3.8. 




sity and dark matter density when the mass in the cells 
exceeded 2x (4x) the mean density in the 0.2 Mpc/h 
(0.5 Mpc/h) simulation. In ah of the GADGET simu- 



■^■^ Our parametrization of the [Morris &: Monaghan| (|l997]> ar- 
tificial viscosity results in a 10 x smaller artificial viscosity with 
ay =0.1 in locations greater than 2 smoothing lengths from where 
this parametrization triggers as having sufficiently larger V • v. 
Our standard implementation assumes ay = I everywh ere, usin^ 
the default GADGET implementation described iniSpringel et al I 

[2MT| 
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and agreement between GADGET and Enzo (see Paper 
II). However, with a fixed gravitational softening length, 
coupling between the dark matter and gas particles spu- 
riously modifies the temperature of the gas as discussed 
in Appendix B. 

We ran the simulations on a combination of comput- 
ing resources including the Henyey computing cluster at 
University of California, Berkeley, and the XSEDE clus- 
ter Trestles. The {0.5, 512^} GADGET simulation re- 
quired 6000cpu-hr on 128 cores of the Henyey cluster 
to reach z ~ 10, whereas the {0.5, 512^} Enzo uni-grid 
simulation required 600cpu-hr to reach z ~ 10. The 
{0.5, 512^} Enzo simulation with four levels of AMR re- 
quired > 10^ cpu-hr on 128 cores to reach z = 20. 

4. RESULTS I: INITIAL FINDINGS 

The top left, top middle, and top right panels in 
Figure [4] show slices of log(l + J^) through three 
{0.2 Mpc//i, 2 X 256^ particle} GADGET simulations 
initialized with v\^^ = 0, v\)c = 30(z/1000) km s~^ [or 
Mach number M^c = 1-9], and ^bc = 60 (z/1000) km s"^ 
[A4bc = 3.8], respectively. Here, is the gas overden- 
sity. Note that 60% of the cosmic volume would have 
had Aihc > 1-9 and 6% would have had A^bc > 3.8 in 
the concordance ACDM cosmological model. The most 
commonly simulated A^bc = case is the least repre- 
sentative: only 10% of space had Mhc < 1- The bot- 
tom panel shows the same slices but through the Enzo 
simulations. Both the GADGET and Enzo simulations 
were initialized with the same random numbers and such 
that the baryons were flowing to the right for the cases 
with Vbc > 0- Note that while the large-scale structures 
largely line-up, there are marked differences between the 
three cases, with larger v\)c leading to a damping of the 
filamentary structures perpendicular to the flow direc- 
tion. The supersonic flows also lead to the development 
of Mach cones around the dark matter substructure. Fig- 
ure [5] illustrates the prominence of these cones, showing 
slices of density in the A4hc = 3.8 simulations. The Mach 
cones are most striking in movies that pan through slices 
in the box. Paper II investigates whether these shocks 
heat the IGM. 

In the remainder of this section, we compare the the 
evolution of the matter in the GADGET and Enzo simu- 
lations both at times where the evolution is linear (such 
that they can also be tested against analytic solutions) 
and also at more advanced stages. 

4.1. Linear Evolution in Simulations 

The most basic test for whether the simulations are 
behaving properly is whether they are capturing linear 
theory at times when it should hold. Linear theory is 
a good approximation for the matter power spectrum in 
the simulations from initialization {z = 200 unless stated 
otherwise) until z = 50, at which point the dark matter 
and baryonic power spectra deviate at the order unity- 
level from its predictions. Figure |6] shows this comparison 
for different initialization methods and codes, comparing 
the power spectrum in the GADGET (left panel) and 
Enzo with AMR (right panel) 0.5 Mpc//i, 2 x 512^ base- 
grid resolution element simulations at z = 150 and 50. 
We find that AMR with our refinement criteria makes 
little difference at these redshifts. The thick curves in 




k [/iMpc-M 

Fig. 6. — Simulated dark matter (top curves) and baryonic power 
spectra (bottom curves) for the cases A^bc = 0? -^bc = l-Q? and 
A^bc = 3.8 (black, blue, and red thick solid curves, respectively). 
The thin solid curves with the corresponding color are the linear 
theory solutions for these A4hc- The curves are calculated from 
the 0.5 Mpc/h, 512^ dark matter particle GADGET and Enzo 
simulations (left and right panel, respectively). The turn-up in 
power at > 10^ Mpc"-*^ owes to shot noise. 

Figure [6] are the dark matter and baryonic power spec- 
tra from these runs, and the thin curves are the pre- 
dictions of linear theory. To calculate these curves, the 
particles in each snapshot were placed onto a Cartesian 
grid with cloud-in-cell interpolation. We then divided 
out the cloud-in-cell window function when calculating 
the power spectrum. The gas density from Enzo was ex- 
tracted onto a fixed grid at its coarsest resolution. We 
then divided out the nearest-grid-point window function 
when calculating each power spectrum. 

At z = 150, both codes trace the linear theory predic- 
tions well, with both the baryon and dark matter compo- 
nents in GADGET and just the dark matter in Enzo de- 
viating at the highest wavenumbers owing to shot noise. 
Whereas, the gas in Enzo underpredicts linear theory 
at wavenumbers that correspond to a few times smaller 
than the Nyquist of the root grid. However, by z = 50 in 
both GADGET and Enzo, there are notable deviations 
from linear theory, especially when A^bc > 0, with the 
dark matter power spectra in Enzo noticeably smaller 
than GADGET at the smallest wavenumbers. 

The thick dashed curve represent the power spectrum 
of the baryons in GADGET for A^bc = when adopting 
the standard practice of running with a fixed gravita- 
tional softening length. Much of this deviation at z = 50 
between the A^bc = linear theory curve and this curve 
is spurious and owes to gas particle-dark matter parti- 
cle coupling. Because the gas overdensity fluctuations 
are just 10~^ at z = 50, even a small amount of nu- 
merical coupling between the dark matter and baryonic 
particles can source significant fluctuations. Appendix 
B provides additional discussion of particle coupling and 
its impact in the simulations. The amount of coupling 
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Fig. 7. — Density power spectra from different simulations at 
z = 20, all with Mhc = 1-9 and using the GADGET code with 
a fixed smoothing length and initialized dX z = 200, unless stated 
otherwise. Top panel: Comparison of our initial conditions gener- 
ator (crosses) with initial conditions using the common practice of 
setting 5i) and 5c to the total matter over density (pluses). Black 
signifies the dark matter power and blue the baryonic power, and 
the curves represent the linear-theory predictions. Second panel: 
The same but instead comparing the simulation with our fiducial 
initial conditions (pluses) with one initialized at 2; = 400 (crosses). 
Note that the two cases largely overlap. Third Panel: Comparison 
of the impact of box size. Shown are our largest particle number 
simulations that are run with box sizes of {0.2, 0.5, 1.0} Mpc//i in 
order of decreasing line width. The dashed curves are the baryonic 
power, and the solid curves are the dark matter power. Bottom 
Panel: Comparison of gas and dark matter power spectra between 
the GADGET and Enzo simulations with AMR, both with specifi- 
cations of 0.5 Mpc//i and 512^ initial gaseous resolution elements. 
The dashed green curve in the bottom panel is the dark matter 
power spectrum from the Enzo simulation, and the dotted red curve 
is the same but instead the baryonic power spectrum. 

does not depend on the box size or particle number of 
the simulation, such that this coupling sneakily evades 
simple convergence tests and can be easily confused with 
nonlinear evolution. This particle coupling is eliminated 
by using adaptive gravitational softening lengths for the 
baryons. This coupling does not only impact simulations 
of the z ~ 20 universe, and, for example, will impact SPH 
simulations of the Lya forest. 

With adaptive softening in GADGET, the deviations 
from linear theory between the GADGET and Enzo 
power spectrum curves largely agree at z = 50. This 
suggests that the simulations in both codes are captur- 
ing linear theory over the realm where it applies and also 
that these deviations are real and linear theory is starting 
to error hy z = 50, especially for A^bc > 0- In addition, 
the matter power spectra in the simulations we ran with 
other box sizes, initialization redshifts, and artificial vis- 
cosity implementations agree well with the simulations 
featured in Figure [6) both at z = 50 and at z = 150. 

4.2. Nonlinear Evolution 

Ultimately, we want to study the nonlinear behavior of 
the gas and dark matter in the cosmological simulations: 



e.g., shocking, collapse of structure, and the formation 
of the first stars. An important step towards this com- 
parison is to verify that the nonlinear solutions of our 
simulations agree. This amounts to comparing the den- 
sity field of the simulations at 2; < 50, times when some 
baryonic density fiuctuations begin to go nonlinear and 
when the first gas-rich halos collapse. Figure [7] shows the 
dark matter (black markers) and baryonic (blue mark- 
ers) power spectra at z = 20. The top panel compares a 
simulation with our initial conditions generator (crosses) 
to one that uses initial conditions with and 5c both 
set equal to the total matter overdensity (pluses), where 
5}) and 5c are the overdensities in gas and dark matter, 
respectively. In addition, the curves in the top panel 
show linear theory, which only captures the growth of 
the larger modes in these simulations. This panel shows 
that the common practice of starting a simulations of the 
first stars with ^5 = 5c results in an overshoot in the size 
of density fiuctuations with a fractional error of Oil). 
This comparison stresses the importance of using differ- 
ent transfer fu nctions for the two components, as also 
emphasized in iNaoz et al. (2011). At higher redshifts 
than are shown, we tind that the power spectrum of the 
case initialized with ^5 = 5c contains significant acous- 
tic oscillations, an artifact from these over-pressurized 
initial conditions. 

The second and third panels down in Figure [7| com- 
pare the convergence between different initialization red- 
shifts and box sizes, respectively. The second panel 
compares a GADGET simulation with our fiducial ini- 
tial conditions initialized at 2: = 200 (pluses) with one 
initialized at 2: = 400 (crosses). Note that the two 
cases are in excellent agreement. However, convergence 
is not as successful in box size as in initialization red- 
shift. The third panel in Figure [7] compares the power 
spectra between GADGET simulations with box sizes of 
{0.2, 0.5, 1.0 }Mpc//i, in order of decreasing linewidth. 
The dashed curves are the baryonic power and the solid 
curves are the dark matter power. The simulations do 
not demonstrate significant convergence with increasing 
box size. We suspect based on calculations discussed in 
Section [2] that if we could run a larger box simulation 
that resolved the Jeans' length, we would conclude that 
the power spectrum in the 1.0 Mpc//i is converged to 
- 10%. 

Lastly, the bottom panel in Figurej6] compares the gas 
and dark matter power between the GADGET and Enzo 
simulations with box sizes of 0.5 Mpc//i and 512^ ini- 
tial gaseous resolution elements. The dark matter power 
spectra largely agree between these two simulations. We 
suspect the small excess in the baryonic power in Enzo 
at /c ~ 10^ Mpc~^ is spurious and owes to power induced 
by the Enzo AMR algorithm. To test this hypothesis, we 
compared our simulations in smaller boxes (with a box 
size of 100//ikpc) with and without AMR. We found 
that for the case without AMR, the power spectra were 
a much better match to the linear theory power spectra 
as the grid size was increased (for z > 50). With AMR, 
we found that spurious extra power induced at the scale 
of the AMR grid migrated to larger scales (smaller k) as 
time proceeded. 

We have focused on the power spectrum as a diagnos- 
tic of the nonlinear evolution, but it is by no means the 
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only measure of structure formation nor the most inter- 
esting one. Subsequently, we will focus on how the first 
structures form in the simulations. While this analysis 
will again provide a comparison of the codes, our focus 
will emphasize new physical insights. 

The most significant improvement of our simulations 
is that the gas is initialized with displacements and 
velocities that are generated from linear solutions that 
include ^bc- This improvement allows us to properly 
follow the infall of gas onto the first cosmic structures. 
The ensuing discussion investigates the impact of v\)c on 
the density distribution of IGM and on the first halos 
that form with deep enough potential wells to retain 
their gas and, in some cases, form stars. The ensuing 
discussion has implications for the thermal evolution of 
the IGM, the formation of the first stars, and the z ^ 20 
21cm signal. 



5. RESULTS II: DYNAMICAL FRICTION IN THE EARLY 
UNIVERSE 

The situation in the early Universe - dark matter 
clumps moving through a more homogeneous gas with 
velocity ^bc ~ is similar to the toy case used to derive 
Chandrasekhar's classic dynamical friction formula of a 
massive particle movi ng through a homogeneous sea of 



collisionless par ticles ( Chandrasekhar 1943 Binney & 



Tremaine 1987|). A similar formula has been shown to 



apply to the case of a par ticle mov ing through gas (e.g., 
Ostrikerl 1999 and Furlanetto fc Loeb,2Q Q4 ) . The dynam- 
ical friction timescale lor a halo of mass to lose its 
energy to the background baryons with density and 
Mach number Mhc 
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where log A is the so called 'Coulomb logarithm'. In the 
gaseous case and for linear perturbations, there exists 
an analytic solution for log A assuming a point-like per- 
turber moving for a finite time ( Ostriker 1999, ): 
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where Rj is the Jeans' scale (or the distance a sound 
wave travels in the age of the Universe f or S = 0, which 
replaces Cg t in the original expression in |Ostriker||1999D . 
To reach equation Q, we replaced rmin (the mini mum 
radius froni which the drag force was included in the|Os 



triker 1999 calculation) with 0.35 A^^-^ times the c. 



riar- 



acteristic scale of an e xtended obje c t. Th is approxima- 



tio n was suggested in Kim & Kim (|2QQ9|) to generalize 



the Ostriker (1999) formula to extended objects sourcing 



nonlinear density perturbations. For definiteness, we set 
this characteristic scale to be Arvir, where rvir(^/i) is the 
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Fig. 8. — Evolution of the component of the baryonic velocity 
in the direction of the baryonic flow times 101/(1 + z). Shown are 
different overdensity cuts in the 0.5 Mpc//i, 2 x 512^ (thick curves) 
and 0.2 Mpc/h, 2 x 512^ (thin curves) GADGET simulations with 
A^bc = 1-9 (-Ubc = 3 km s"-*^ at z = 100). The top panel shows 
the average velocity above the specifled overdensity thresholds at 
the specifled redshift. The bottom panel follows the velocity of a 
group of gas particles that fall with in the same density range at 
the time of initialization. There are 52, 15, 21, 9, and 3% of the 
particles in each density grouping specifled in the bottom panel for 
the 0.5 Mpc/h box. Both panels illustrate how dynamical friction 
causes the gas and dark matter to decelerate into the same frame 
in overdense regions. 

virial radius of a halo of mass mh- Note that with equa- 
tions Jsl) and (Ig]), tdf is at a minimum for A4hc = 1 — 2 
(OstriEer^^l999 j7 which coincides with the most probable 
values lor Mhc in the concordance cosmological model. 

A Hubble time equals 280 Myr at times when 

Qrn ^ 1- Thus, equation ^ implies that the surrounding 
gas will be decelerated and fall into halos with nih ^ 
10^~^ Mq by z ~ 20. This bound coincides roughly with 
the mass bound on halos that can ove rcome pressure and 
maintain their gas ( Naoz et "ar]|2009 ). 

Significant dynarnical Iriction occurs when nonlinear 
structure forms such that gravity becomes strong enough 
to generate 0{1) perturbations in the gas on smaller 
scales than the distance a sound wave travels in a 
Hubble time. The fact that log A ~ 1 refiects that there 
is not a large range of scales between the sizes of objects 
and how far sound waves can travel. 

Gaseous dynamical friction generates Mach cones 
around supersonic objects, in our case predominantly 
dark matter halos. The shocks at the edges of these 
cones may heat the IGM as the halo is decelerated. These 
cones are very visually apparent in the large-scale sim- 
ulations with finite Vhc (see Fig. [s]) and especially when 
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one zooms in on halos (as we will show in Figure [s] 
plots the evolution of the average velocity of gas particles 
above a fixed overdensity threshold (top panel) and the 
average velocity of particles that fall into the same den- 
sity range at initialization (bottom panel) in GADGET 
simulations with A^bc = 1-9. Both panels show that dy- 
namical friction effectively acts to decelerate the densest 
regions. Interestingly, the velocity even goes negative for 
the rarest, most overdense regions in the bottom panel, 
as the wake of baryonic material around a dark matter 
potential halts and falls back onto the dark matter poten- 
tial well. We find very similar trends in the simulations 
with Aihc = 0.6 and 3.8 as in the Aihc = 1.9 simulations. 

This dynamical fictional also causes the velocity of 
the dark matter halos to move in the direction oppo- 
site the decelerating baryons. However, because there is 
^ 8 times more mass in the dark matter, the change in 
velocity of dark matter overdensities from this dynamical 
friction is <C ^bc- 

6. RESULTS III: THE IMPACT OF Vbc ON THE FIRST 
STRUCTURES AND HALOS 

A supersonic dark matter-baryon velocity difference 
affects how much gas falls onto each halo as well as 
the profile of intra-halo gas. For relatively massive ha- 
los (> 10^ ^o), the density of the gas in turn impacts 
whether they can form enough molecular hydrogen to 
cool quickly enough and form stars. However, it is not 
clear that such massive halos should be significantly im- 
pacted by Vbc- Halos with > 10^ Mq have circular ve- 
locities of Vcir > 2.5 km s~^ at z = 20, whereas the 
RMS value for ^bc is 0.6 km s~^ at that redshift. We 
investigate this question here. We show that the impact 
of v\)c is actually quite significant, even for halos with 
Vcir = 7 km s~^ (corresponding to masses of 3 x 10^ ^o)- 
Analytic arguments for why this should be the case are 
presented in Paper II. 

Figures [9}|TT] compare density slices of three individ- 
ual halos at z = 20 in simulations with {0.5Mpc//i, 
512^ initial gaseous resolution elements}. The top row 
of panels use GADGET and the middle Enzo.^^ The 
bottom row of panels shows the density contrast of the 
baryons to the dark matter or (1 + ^5) /(I + ^c)- Fig- 
ures [9] and [To] show the two largest halos in these sim- 
ulations at z = 20. Their friends-of-friends group 
mass in the GADGET simulation with Vbc = are 
2 X 10^ Mq and 8 x 10^ Mq (using a linking length 
of 0.2). These masses are slightly smaller in the simula- 
tions with nonzero ^'bc- Figure pT] zooms in on a slightly 
lower mass halo, the ~ 20*^ most massive halo in the 
simulation {rrih ~ 2 x 10^ Mq). This halo was selected 
(out of hundreds of other halos) because its environment 
demonstrates a large variety of effects owing to the rela- 
tive velocity. 

The visual morphology of the virialized region in the 
most massive halo in the box is not changed significantly 
between the A4hc = and A^bc = 1.9 cases, but the 
morphology of this halo for A^bc = 3.8 is significantly 
changed (Fig. [9|. Quantitatively, the central density is 
decreased by a factor of 2 between the Mhc = and 

12 The GADGET slices are generated using th e SPH kernel. The 
Enzo slices are generated using the yt package ( [Turk et al.pOll]) 
and such that the cell thicknesses depends on Ihe level ot AlVLK 
refinement. 



A4hc = 1-9 cases in Enzo, and by more than a factor of 
5 between the Mhc = and Mhc = 3.8 cases. In the 
GADGET simulations, the magnitude of this suppres- 
sion is similar to that in Enzo, but the central density is 
always smaller in the GADGET simulations. The differ- 
ences in central density becomes somewhat larger as the 
halo mass is lowered (see 



6.2) 



In addition, the accretion of gas is impacted by ^bc for 
the concentrated filament that forms perpendicular to 
the bulk fiow in Figure 9] The orientation with respect 
to i;bc of fiows onto halos has a notable impact on the 
amount of gas accreted by the halos in our simulations, 
and these effects introduce significant stochasticity in the 
baryonic mass fraction of halos at a fixed dark matter 
mass. The bottom panel in Figure |9] shows that the dark 
matter is still present in these filamentary fiows, but the 
baryons are gone, with a slight overdensity downwind. 

The second most massive halo (Fig. 10) is even more 
altered by the bulk fiow of the baryonsior Mhc = 1.9 
than the most massive halo. (Note that there is a neigh- 
boring halo that appears in these panels with a somewhat 
smaller dark matter mass.) Most striking is the appar- 
ent bow shock that develops near the virial radius of this 
halo (and its neighbor) as well as the wisps of down- 
wind gas. Remarkably, the halo has little bound gas in 
the case with Mhc = 3.8. Again, we see that filamentary 
structure perpendicular to the baryonic flow is disrupted. 

Lastly, Figure pT] shows a halo that is an order of mag- 
nitude less massive, with rrih ^ 10^ Mq. In this case, 
Mach cones develop around many of the dark matter 
overdensities and not just the central halo. While this 
halo is just at the mas s threshold that is cap able of cool- 
ing and forming stars ( Machacek et al.|2001 ), there is no 
way it can form a star in the Mhc = 3.8 case as the den- 
sity of the contained gas is only an order of magnitude 
above the cosmic mean. 

For the case with Mhc = 0, the baryons and dark 
matter largely trace each other in the simulations in pro- 
portion to their cosmological abundances on scales larger 
than the Jeans' scale of the gas. When there is a nonzero 
velocity difference between the dark matter and baryons, 
the dark matter and baryons no longer trace each other 
even in the most massive halos in our simulations (see the 
bottom panels in Figs. ^11). Even the highest density 
gas in the most massivenalo is offset from the center of 
the dark matter halo when Mhc = 3.8.^^ In general, the 
gas trails downwind of the dark matter, especially when 
the filamentary structure is aligned perpendicular to the 
bulk fiow of the gas. Indeed, it is in these filaments where 
the gas is often completely swept out of the dark matter 
well. The offset of baryonic mass delays accretion as the 
baryonic wake turns around and falls back onto the halo. 
We find that its fallback velocity generally exceeds the 
typical velocity at which gas is accreted onto a halo in 
the case Mhc = 0, and the fallback of course is lopsided, 
which further buffets the halo and reshapes the internal 
gas distribution. We also found that nonzero ^bc can re- 
sult in the total angular momentum of the gas changing 
direction. 

In Figure [12] investigates how ^bc impacts the forma- 

We note that because we have centered the images on the 
densest peaks, each frame tends to move with the baryon fiow as 
A^bc increases. 
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Fig. 9. — Comparison of the gas density in the most massive halo in the {512^, 0.5 Mpc/h} simulations at z = 20. Plotted in the top 
(middle) row is a cut through this halo with width 50 kpc/h in the GADGET (Enzo) simulation. From left to right, A4hc = 0, 1.9, and 
3.8, with the baryons moving to the right. This halo's friends-of-friends group mass in the A^bc = GADGET simulation is 2 x 10^ Mq. 
Plotted in the bottom row is the density contrast of the baryons to the dark matter or (1 + + Sc) in the GADGET simulation. Each 

slice is centered on the densest point within the halo. 



tion of shocks around the two most massive halos in our 
simulations at z = 20. The shocks that heat the gas are 
largely confined to the immediate region around the ha- 
los, especially when ^bc > 0? and only the strongest Mach 
cones generate appreciable entropy. At z = 25, when 
these halos are beginning to virialize, the A4hc = 1.8 ha- 
los have bow shocks even before the virial shocks develop 
in the A4hc = 0.0 halos. 

We have also analyzed the impact Vbc has on the vortic- 
ity of halos gas. Vorticity is defined as V x v^, where is 
the gas velocity field. The cosmic initial conditions have 
zero vorticity, but vorticity is required to seed turbu- 
lence and magnetic fields (see e.g., Pudritz fc Silk^l989 ), 
which may alter the properties of the first stars, 'l^e pri- 
mary way to generate cosmological vorticity is via curved 
shocked fronts, and ^bc clearly alters the morphology of 



shocking regions. We found that there was only a mod- 
est difference in the magnitude of voriticity with 'Ubc > 
in the three most massive halos at z = 20, when most 
of the gas had already decelerated. However, at z = 25 
the most massive halos in simulations with A^bc = 1-8 
had an order of magnitude more vorticity compared to 
those in that with Mhc = 0.0. At z = 20, the - 10^ 
halos in the simulations had higher vorticities than in 
the case with 'Ubc > 0. (Note that while cosmic vortic- 
ity should be conserved, it is almost certainly damped 
out with time by the high artificial/numeric viscosity in 
our simulations.) These trends suggest that for 'Ubc > 0, 
halos develop shocks earlier and have greater vorticity 
(and, thus, are more turbulent) than halos in previous 
simulations of the first halos and stars with v\yc = 0. 
Overall, the GADGET and Enzo simulations show a 
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remarkable level of agreement on the structure and mor- 
phology of the three halos in Figures [9- 11 The location 
and shape of the gas and dark matter structures are sim- 
ilar in both codes' simulations on scales as large as the 
box size 1 Mpc; Fig. |4| to as small as a tenth of the 
virial radius of 10^ M© halos (- 0.1 kpc; Fig[l0|. The 
halos in GADGET tend to be less dense than tHe halos 
in Enzo. One possible cause of this may be the different 
effective resolutions between two codes: The minimum 
softening length in GADGET at z = 20 is 0.04 kpc (ap- 
proximately two decades smaller than the halo diameter) , 
and the AMR refinement criteria that was used for Enzo 
is super-Lagrangian. However, in testing Enzo simula- 
tions with smaller boxes and more aggressive refinement, 
we found that the earlier the refinement was initialized, 
the more concentrated the halos appeared at later times. 
Note that none of the simulations includes molecular hy- 
drogen cooling and so this comparison is not sensitive to 



cooling rates. The simulations also have similar ther- 
mal properties throughout the volume, again with the 
exception of the peak densi ties f ound in the halos. We 
refer the reader to Appendix |B.2| to directly compare the 
thermal propertiers of the gas in the various simulations, 
and the impact particle coupling has on the GADGET 
simulations. 



6.1. The First Structures 

As is evident from Figures [9) [lO) and [TT] the infall of 
gas onto the first halos is suppressed in the A^bc > 
simulations relative to the Alhc = _Q runs for typical 
streaming velocities ( Stacy et al | 2 01 1 [ |Greif et al.|2011b 



Fialkov et al. 2011). In 1^'igure |13[ we investigate the 
amount of mass in baryons that accretes onto each halo 
as a function of the dark matter halo mass for Enzo 
simulations with 0.5 Mpc//i, 512^ initial gaseous reso- 
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Fig. 11. Same as Fig. [o] 
GADGET simulation is 2 x 1 



but shown is the 
0^ Mq. 



20^^ most massive halo at z = 20. Its friends-of-friends group mass in the A^bc 



lution elements. Below 10^ Mq in this plot, we instead 
show the logarithmic mean of the baryonic mass and its 
variance for the halos within a mass bin. (The phys- 
ical mean of these curves is similar to the logarithmic 
mean.) There is no well defined way to measure the bary- 
onic mass contained within the first structures, especially 
when A^bc > so that the dark matter and baryons of- 
ten do not trace one another. In many instances (even 
for the most massive halos in our simulations), the gas 
and dark matter density peaks are offset by as much as 
r2oo - the distance from the halo density peak at which 
the the dark matter overdensity falls below 200 - so that 
it is unclear what is the best method to center when re- 

The dark matter mass can be decreased by as much as a factor 
of two for the case Ai^c = 3.8, even in the largest halos. In Paper 
II, we discuss in detail the suppression of the dark matter mass 
function for the largest GADGET simulations when A4bc = and 
1.9. 



porting spherically averaged quantities. We choose to 
center on the gas density peaks in this study. To find 
these peaks, we first search for the densest peak within 
two viral radii of the dark matter center of mass, and 
then we select a sphere with radius r/^, such that the 
total enclosed dark matter mass is rrih, as found by the 
HOP halo finder dEisenstein fc Hut|[l998| ). While this 
algorithm is non-traditional, it yields a conservative es- 
timate for the amount of mass that is contained. For 
example, we find that if we use a criterion more simi- 
lar to that used in other analyses such as Na oz et ciL] 
(2011) (in which we start with the dark matter density 
peak or dark matter center of mass, or especially if we 
calculate the gas within r2oo rather than r^), the result- 
ing baryonic masses are significantly lower for halos in 
simulations with A4hc > 0- In addition, the halo-to-halo 
scatter also increases. For example, if we center on the 
dark matter center of mass, as opposed to gas density 
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Fig. 12. — Impact of i;bc on shocks in the two most massive halos. Plotted are shces of the temperature a nd t he "entropy," which scales 
as 

T/p2/3^ of the most massive halo (top) and second most massive halo (bottom) as shown in Figs. [9|and |l0| respectively. The entropy 
generated from the shocks is confined to the immediate vicinity of the halo. Note: The frame size shown is a tact or of two smaller than in 
Figs. |9] and [To] to show more detail. 
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Fig. 13. — Mass of baryons in halos as a function of halo mass, 
ruh- For each halo, we find the peak gas density and, then, de- 
termine the mass of baryons enclosed within r^, the radius that 
encloses of dark matter. We describe the reasons behind this 
nontraditional setup in the text. The black circles, red stars, 
and blue triangles represent halos from Enzo simulations with 
0.5 Mpc/h, 512^ initial gaseous resolution elements and for the 
cases A4hc = O-O? l-Q? and 3.8, respectively. For ruh > 10^ Mq, we 
show each individual halo. The remaining data at < 10^ Mq 
show the mean and standard deviation of the logarithm of the val- 
ues for individual halos. The dashed line equals Qb/Qcmh- 

peak, and calculate the gas within a distance r/i, this de- 
creases the baryon mass by 0.94 ± 0.16, 0.73 ± 0.38, and 
0.50 ± 0.36 on average (zb the standard deviation of the 
scatter) for A^bc = O-0, 1-9, and 3.8, respectively. Fi- 
nally, we find consistent results between the Enzo and 
GADGET simulations using a similar gas-finding algo- 
rithm. 

Figure [13] illustrates that increasing the value of A4hc 
in the simulations significantly reduces the mass of 
baryons that have accumulated into dark matter ha- 
los. To a lesser extent, it also lowers the mass of dark 
matter that has collapsed into the most massive halos. 
In addition, the scatter in the baryonic mass at fixed 
halo mass increases drastically with A^bc- For a 10^ Mq 
halo and A4hc = 3.8, the baryon fraction is approxi- 
mately half that found when there is no relatively ve- 
locity. The halo-to-halo scatter in the baryonic fraction, 
on the other hand, increases to 100 % from only 10 % be- 
tween At be = and A^bc = 3.8. For cases in which there 
is a bulk velocity between the dark matter and baryons, 
the enclosed baryon mass in individual halos becomes 
dependent on the environment of the halo. From our 
earlier investigations of individual halos, it is clear that 
halos found in filaments perpendicular to the bulk fiow 
appear more affected, as do halos without any apparent 
nearby structure. Such environmental dependences are 
likely sourcing the added stochasticity as I'bc is increased. 

6.2. Star Formation 



Fig. 14. — Mass of baryons in halos that can cool in a Hubble 
time as a function of halo mass in the 0.5 Mpc/ Enzo AMR simu- 
lations. For each halo, we find the mass of baryons enclosed in 
that can cool in a Hubble time, where is the radius that encloses 
rrih of dark matter. The black circles, red stars, and blue triangles 
represent data from the A4hc = O-O? l-Q? and 3.8 respectively. The 
dashed line equals Qb/^cTrih- 

As we described in § |6.1[ the baryon-dark matter veloc- 
ity difference suppresses the accretion of gas into the first 
halos, thereby reducing the density and temperature of 
these structures. This suppressed accretion also strongly 
damps the rate of star formation in the first halos. Our 
simulations do not follow the gas to densities where it 
becomes a star since we do not have molecular hydrogen 
cooling. Thus, we use the mass of gas that can cool in a 
Hubble time from molecular hydrogen cooling as a proxy 
for the star formation rate in our simulations. However, 
we find that our relative comparisons among the simula- 
tions do not change if we instead take the gas mass that 
can form within 0.1 Hubble times. In our opinion, it is 
not necessarily a disadvantage to use our cooling crite- 
ria as a proxy for star formation rather than following 
the cooling and condensing gas to much higher densities: 
Firstly, simulations that follow gas to much higher den- 
sities do not all agree o n the character of star formation 



in the first halos (e.g., Greif et al.| 2011a). In addition 



the halo (Alvarez et al.||2006 


Yoshida et al.|2007D or the 


cosmological Lyman- Werner 


background (Haiman et al.| 



plexity of modeling such star formation beyond when the 
first star has formed in the halo or when the cosmolog- 
ical star formation rate density has reached a value of 
- 10-^ M^yr Mpc"^ 



Figure 14 shows the mass in baryons that have cool- 
ing timesless than a Hubble time as a function of 
the halo mass for our Enzo simulations with boxsize 
0.5Mpc//i and 512^ initial gaseous resolution elements. 
To calculate th e cooli ng time, we use the formula in 



Tegmark et al. (1997) under the crude approximation 
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Fig. 15. — Peak gas density in halos as a function of mass in the 
0.5 Mpc//i, 512^ base grid Enzo AMR simulations. For each halo, 
we find the peak gas density within two viral radii of the halo center 
of mass. The black circles, red stars, and blue triangles represent 
data from the A^bc = O-O? l-Q? and 3.8 simulations, respectively. 
For rrih > 10^ Mq, we show the peak density for each individual 
halo. The remaining points with error bars shows the mean and 
standard deviation of the logarithm of the values. For a given dark 
matter halo mass, the peak density is reduced by approximately 
half 1/8) when A^bc = 1-9 (3.8) as compared to Mhc = 0.0. 
At z = 20, the mean gas density is 4 x 10~^^ g/cm"*^. 

that log(l + A/'rec) = 1 to calculate the amount of molec- 
ular hydrogen, where A/'rec is the number of recombina- 
tions. This equality approximately holds for gas at the 
virial density of halos at z = 20, and it allows us to 
calculate the cooling rate from a single time slice rather 
than following the density evolution of a fluid element. 
To tabulate the amount of cooling gas, we search within 
a sphere of from the halo center of mass. However, 
we find that our qualitative results are insensitive to this 
choice and searching within r2oo around the gas density 
peak yields similar numbers. 

We find that the dark matter mass threshold required 
to have high enough temperatures and densities to cool 
increases with A^bc in agreement with the picture in pre- 
vious studies (|Greif et al. 2011b; Stacy et al. 2011J [Fi-] 
alkov et ah 201 ip . Surprisingly, we find that even some 
of the most massive halos in our boxes (up to 10^ Mq in 
the 1 Mpc boxes reported in Paper II) have less cooling 
gas for the case with Mhc = 1-9 than with Mhc = 0, 
in agreement with |Naoz et al.| |2011| . The amount by 
which the star formation rate is decreased has signifi- 
cant variance at a single mass. On a halo-by-halo basis 
{98%, 57%, 11%} of halos with dark matter halo masses 
rrih > 4x 10^ M©, have more than 100 Mq of gas that will 
cool in a Hubble time for Mhc = {0.0, 1.9, 3.8}. Even 
for ruh > I X 10^ Mq, ^ 2/3 of halos in the Mhc = 3.8 
simulation have less than 10^ Mq in gas that can cool. 
We just note here that Figure p!4f uses our Enzo simula- 
tions, with the HOP halo finder to determine dark matter 



halo sizes, which we found to be reduced compared to the 
dark matter halos in GADGET. 

The impact Mhc has on the peak gas density in ha- 
los is related to the impact Mhc has on the amount of 
gas that can cool. In Figure [Tsj we show the peak gas 
density in each halo in the 0.5 Mpc//i, 512^ base grid 
Enzo AMR simulations as a function of the halo dark 
matter halo mass. We find that the peak central density 
decreases by nearly an order of magnitude, on average, 
when Mhc = 3.8. The dispersion of peak densities also 
increases significantly, similar to what we found in the 
baryon masses of the halos. 

In Paper II, we exclusively use the GADGET simula- 
tion on a larger box to look at how the total star for- 
mation rate is modulated by Mhc In the smaller box 
sizes reported here, we find that Enzo, on average, has 
higher gas densities and temperatures in the halos (see 
our discussion in SS 



B.2 



and Fig. |16[ ) . This results in ap 
proximately four times as much gas that can cool to form 
stars. In addition, the number of halos that form stars 
is approximately twice as large as found in GADGET. 
However, the fraction that star formation is reduced by 
Mhc is the same in both codes. 

As a final aside, given the large offset between the 
baryons and dark matter in our simulations, it seems 
improbable that the first stars ever form at the dark mat- 
ter density peak in the halo. Thus, we find it u nlikely 
that stars supported by dark matter annihilation (|Spol- 
|yar et al.H 2008) ever existed, even with favorable dark 
matter annihilation cross sections. 

7. CONCLUSIONS 

This study investigated how gas was captured into the 
first structures during the cosmic Dark Ages. This cap- 
ture is relevant to the formation of the first stars and the 
high-redshift 21cm signal. In the concordance ACDM 
cosmological model with nearly scale-invariant primor- 
dial potential perturbations, the initial conditions that 
set the distribution and velocities of matter are well con- 
strained by observations of the CMB and of large-scale 
structure. This study aimed to follow the evolution of 
these (initially small) matter fiuctuations into their non- 
linear state at lower redshifts. However, there are chal- 
lenges with achieving this objective. A simulation must 
have the resolution to capture scales below the Jeans' 
scale of the gas as well as be large enough to capture a 
representative sample of the Universe. We discussed the 
specifications necessary to meet these competing objec- 
tives. 

We ran a large number of simulations of the early 
Universe with both the GADGET and Enzo cosmology 
codes. Our simulations are the first to initialize the 
gas and dark matter self-consistently from linear theory 
on scales where gas pressure is important and at red- 
shifts where the baryons and dark matter do not trace 
each other. There is no physical or significant compu- 
tational reason to not initialize cosmological gas+dark 
matter simulations with anything but the full linear solu- 
tions, but this appears to not have been done previously. 
In contrast to previous numerical studies, we initialized 
most of our simulations in a manner that included (1) 
the impact of pressure on the growth and rate of growth 
of modes, (2) temperature fiuctuations, (3) the correct 
baryonic and dark matter velocities, and (4) transfer 
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functions that account for the dark matter-baryon veloc- 
ity differential (under the approximation that Compton 
drag from the CMB was zero at z < 1000). We showed 
that simply boosting ^bc at the initial redshift of the 
simulation, as done in prior studies, misses much of the 
impact of v\)c on the linear growth of density fluctuations. 

Because our simulations were run with two very 
different codes (the tree-particle mesh plus SPH code 
GADGET and the nested particle mesh plus AMR code 
Enzo), the comparison of these codes' simulations pro- 
vided a relatively controlled test of the codes themselves. 
Despite their very different algorithms, the ~ 10^ Mq 
halos that formed in GADGET and Enzo (with AMR) 
simulations have similar appearances and their thermal 
properties agree remarkably well. While the densities 
in these halos were lower in GADGET, we found no 
evidence that GADGET either significantly under- or 
over-estimates the entropy produced in weak shocks 
from structure formation despite its SPH hydro-solver. 
The largest difference between the simulations with the 
two codes was from spurious gas-dark matter particle 
coupling in the GADGET simulations, which was most 
pronounced when ^bc = 0- We found that this coupling 
is slightly smaller when we use the common prescription 
of staggering gas and dark matter particles at half an 
interparticle spacing (but with the same glass file) as 
opposed to using two separate glass files for the gas and 
dark matter. This coupling is likely to have impacted 
all previous cosmological SPH simulations. However, by 
using adaptive gravitational softening for the baryons 
in the simulations, we were able to eliminate the bulk 
of this particle coupling. We recommend that all 
cosmological SPH simulations use adaptive gravitational 
smoothing lengths, as the effects of particle coupling 
evade standard convergence tests. 

Our primary focus was to model the impact of v\)c dur- 
ing the cosmic Dark Ages on the formation of the first 
nonlinear structures, especially the first stars. To do so, 
we used a suite of ~ 20 GADGET and Enzo simulations 
to explore the formation of the first structures in the 
Universe as well as the thermal evolution of the gas. We 
found that: 

1. The halo gas in these simulations is often found 
significantly downwind and with lower densities in 
the simulations with Mhc > 0- The lower densities 
delay the formation of the first stars. This delay is 
consistent with what was found in prev ious studies 



(Stacy et al. 2011 Greif et al. 2011b). These ef- 



lects impact the lormation ot the first stars and 
could imprint fluctuations in high-redshift 21cm 
backgrounds (as discussed in Paper II). The gas 
that accumulates in the halos is off-center from the 
dark matter halo, with the density peaks in the 
gas and dark matter offset by as much as r2oo foi" 



the case A^bc = 2 (which is near the cosmological 
average of A^bc)- Furthermore, we find that the 
maximum gas density can be suppressed on aver- 
age by roughly an order of magnitude between the 
cases A^bc = and A^bc = 4 for halos with masses 
of 10"^ Mq and even up to 10^ Mq. 

2. Dynamical friction induced by the gas flowing by 
the dark matter halos induces visible Mach cones, 
and the tug from the mass in these supersonic 
wakes acts to erase the velocity difference between 
the dark matter and baryons. We showed that the 
dynamical friction timescale for > 10^ Mq halos 
is less than a Hubble time for typical At be- Most 
of the overdense gas in our simulation has deceler- 
ated into the dark matter frame via this process by 
z = 20. These downwind wakes eventually crash 
back onto dark matter halos, with larger infall ve- 
locities than most of the accreted gas. However, we 
found that these shocks did not significantly heat 
the IGM. 

3. The halo environment plays a significant role in 
how much gas makes it into a halo in cases with 
typical Albc- Eoi" example, we find that dark mat- 
ter filaments perpendicular to the gas flow are often 
devoid of gas. Whether a halo has gas or not de- 
pends on the orientation of intersecting filaments. 
As a result, the baryonic mass fraction of halos (as 
well as the fraction of mass that can cool and form 
stars) shows significant stochasticity at fixed halo 
mass in the simulations with Albc > 0- 
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APPENDIX 
A. LINEAR THEORY 

The system of equations for the hnear evolution of the matter in the presence of an initiahy uniform velocity difference 
between the dark matter and baryons is ( [Tseliakhovich fc Hirata||2QlQ| ) 

^c = -^c, (Al) 

Oc = -3H^/2 {n^Sc + ntSb) - 2H0c, (A2) 

6b = -i Vbc ■ kSb- Ob, (A3) 

Ob = -i a"^ Vbc -kOc- 3H^/2 {njc + ^bSb) - 2H0b + k^a-'^Sb, (A4) 



where Ox is related to the velocity field via Vx = —iak ^kOx (since vx has zero vorticity until the onset of non- 
linearity for infiationary initial conditions), and c denotes cold dark matter and b baryonic m ater ia l. T hese equations 



use that the homogeneou s component of vpc r edshifts away as 'Ubc oc (1 + 2;). Equations (Al - A4) use a similar 
notation to that in Tselia khovich fc Hirata|[2QTQ except that they are in the frame of the dark matter (which results 
in less oscillatory solutions m imagmary space once the baryons begin to fall into the dark matter potenti al wells). 
An addit i onal e quation for the temperature is required to solve for Cg. We use the temperature equation in |Naoz fc 



Barkana (2007), which includes adiabatic processes and Compton heating off of the CMB. 

B. PARTICLE COUPLING 

During the initial comparison of the GADGET and Enzo simulations, we found a significant excess in power in 
the gas density power spectrum on small scales in our GADGET simulations relative to Enzo (and to linear theory), 
especially in the GADGET simulations with '^bc = (see Fig. [6|. This exces s pow er remained with the same mag- 
nitude and slope even when the particle resolution was increased. Appendix |B.1| describes how this excess owes to 
particle coupling, and Appendix |B.2| shows how the coupling alters the distribution of temperatures and densities 
in our GADGET simulations. For the simulations discussed in the body of this paper, we used a variable gravita- 
tional softening length (with this set equal to the softening length), which successfully removes at least the bulk of 
this coupling. Variable softening is not typical in such simulations, but should be done in all further cosmological 
SPH simulations. It simply involves turning on the compile-time fiags ADAPTIVE_GRAVSOFT_FORGAS (and also 
ADAPTIVE_GRAVSOFT_FORGAS_HSML in GADGETS). 
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B. 1 Origin of Coupling 

Here we attempt to understand how particle coupling should scale with redshift, box size, and particle number. This 
boils down to comparing the escape velocity of particles from each other's potential well to their relative velocity. The 
escape velocity at a proper distance r from a dark matter particle of mass Mp is 

?ii^4V" = o^63('^V'V^)"'feV"' km.-', (Bl) 



bM^J V2007 

where /ms is r in units of mean interparticle spacings and note that 5 Mq is Mp in the 0.2 Mpc//i, 2 x 512^ particle 
GADGET simulations. The gravitational force increases down to the softening length, which we have taken to be 

/ms = 0-04 in our GADGET simulations. (The rule of thumb for cosmological simulations is 0.03 0.04.) 

Compare Vesc,p to the relative velocity of particles. The Hubble flow dominates over peculiar velocities in determining 
the relative velocity of neighboring particles. The Hubble velocity of particles during matter domination is 

vn = Hr = ^mi(^\^' (^]'^"(^\ kms-. (B2) 

V200/ \omJ ^ ^ 

Both Vesc,p and vh have the same scaling with redshift and Mp. Thus, box size and particle number do not change 
the typical amount of coupling between a particle and its neighbors, explaining the trends seen in our numerical tests. 
Also, Vesc,p and vh are only comparable at approximately one interparticle separation. At smaller separations, Vesc,p 
is larger, meaning that particles are likely to become trapped in the potential well of other particles. 

However, the velocity of the baryons relative to the dark matter naturally alleviates some of this coupling as 
'^bc ~ 5 X (z/200) km s~^, so that the energy of gas particles is such that they can travel in and out of the potential 
wells of the dark matter particles. However even in the absence of dynamical friction, the different scaling with z 
between the relative velocity and the escape velocity ensures that even in the case with ^bc > 0, some particle coupling 
is unavoidable with a fixed smoothing length. 

B.2 Impact on Temperatures and Densities in the Simulations 



In Figure 16, we plot temperature-density phase diagrams for six simulations. The top (bottom) panel shows a 
two-dimensional histograph of gas mass as a function of temperature and density for the GADGET (Enzo) simulations 
with A4hc = 0.0,1.9, and 3.8 (from left to right, respectively). The top right panel shows the impact of part icle 



coupling. The region above the dashed lines can cool in a Hubble time owing to molecular hydrogen (see ^6.2 for 
details), and the dotted lines show a single arbitrary adiabat {Tg ex p^^^). The deviation from a single adiabat at 
low temperatures and densities is determined by Compton heating off of the cosmic microwave background. Although 
subtle, near the mean density of the gas (p ^ 4 x 10~^^g/cm^) the simulations with a larger dark matter-baryon 
velocity difference have an ~ 1.5 times broader distribution in temperature at fixed density. 

In Figure[T7| we show probability distribution functions (PDFs) of the gas density and temperature in the simulations. 
For v\yc > OTTne gas density PDFs in both the Enzo and GADGET simulations largely agree, with some gas in the Enzo 
simulations reaching higher densities than in GADGET, likely because Enzo has higher resolution in over densities. 
For v\)c = 0, we find significant disagreement between the overall shape of the density PDF, as highlighted in pink in 
Figure pT} This disagreement owes to spurious particle coupling between the dark matter and gas in GADGET. The 
amount of particle coupling is somewhat reduced when staggering the gas particles at half an interparticle spacing in 
each dimension from the dark matter particles (whose locations are set by a g lass file as is often do ne; blue dash-dotted 



curves) rather than using two separate glass files (the method advocated in |Yoshida et al.||2003[ blue dotted curves). 
This reduced coupling results because the dark matter and gas particles have a reduced probability of being near one 
another in the staggered case. Thus, we also recommend staggering particles rather than two separate glass files. 

The impact of particle coupling is especially apparent when comparing the peak densities that the GADGET particles 
reach in the different simulations. When ^bc > 0? the GADGET simulations appear to be less dense than Enzo at the 
highest densities. For v\yc = 0, the simulations with GADGET have denser gas than the Enzo simulations. There are 
also significant disparities at ^5 ~ 10 {pb ^ 10~^^ g cm~^) in the 'Ubc = case between the GADGET and Enzo. 
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Fig. 16. — Temperature-density phase diagram of the z = 20 gas in the 0.5 Mpc//i, 2 x 512^ base-grid resolution element simulations. 
The top (bottom) panel shows a two-dimensional histograph of gas mass as a function of temperature and density for the GADGET 
(Enzo) simulations with A^bc = 0.0, 1.9, and 3.8 (from left to right, respectively). The top right panel should be compared to the top left 
panel. It shows the effect of particle coupling on the phase diagram from the GADGET simulation that used a fixed softening length with 
-^bc = O-O- For reference, the region above the dashed lines can cool in a Hubble time owing to molecular hydrogen (see the text for 
details), and the dotted lines show a single arbitrarily normalized adiabat (temperature oc p^/^). Each point is color coded to represent 
the mass that falls in one of 400 x 400 bins in log space. GADGET contains none of the lowest mass, purple points because it has a gas 
particle mass of 16.75 M©. 
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Fig. 17. — Histograms of the density (left panel) and temperature (right panel). In both panels, we plot mass-weighted PDFs with solid 
(dashed) curves from our 0.5 Mpc/h GADGET (Enzo) simulations. The {black, red, green} curves correspond to A^bc = {0.0,1.9,3.8} 
(roughly from top to bottom). The density histograms for Enzo and GADGET agree rather well at < 10~^^g/cm^. The blue dotted 
curves show the results of our simulations with A^bc = 0-0 cind a fixed gravitational softening length. We attribute the discrepancy between 
these curves and the other GADGET curves (solid black) to gas-dark matter particle coupling. The region where this disparity is greatest 
is highlighted in pink. The blue dash-dotted curves show the results of the Albc = 0-0 simulations where gas particles are staggered at half 
a mean particle spacing in each dimension from the dark matter particles. 



